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Abstract. Stochastic spatial predator-prey competition models represent paradigmatic systems to under- 
stand the emergence of biodiversity and the stability of ecosystems. We aim to clarify the relationship 
and connections between interacting three-species models and the classic two-species Lotka-Volterra (LV) 
model that entails predator-prey coexistence with long-lived population oscillations. To this end, we utilize 
mean-field theory and Monte Carlo simulations on two-dimensional square lattices to explore the tempo- 
ral evolution characteristics of two different interacting three-species predator-prey systems, namely: (1) 
a cyclic rock-paper-scissors (RPS) model with conserved total particle number but strongly asymmetric 
reaction rates that lets the system evolve towards one "corner" of configuration space; (2) a hierarchical 
"food chain" where an additional intermediate species is inserted between the predator and prey in the 
LV model. For the asymmetric cyclic model variant (1), we demonstrate that the evolutionary properties 
of both minority species in the (quasi-)steady state of this stochastic spatial three-species "corner" RPS 
model are well approximated by the two-species LV system, with its emerging characteristic features of 
localized population clustering, persistent oscillatory dynamics, correlated spatio-temporal patterns, and 
fitness enhancement through quenched spatial disorder in the predation rates. In contrast, we could not 
identify any regime where the hierarchical three-species model (2) would reduce to the two-species LV 
system. In the presence of pair exchange processes, the system remains essentially well-mixed, and we 
generally find the Monte Carlo simulation results for the spatially extended hierarchical model (2) to 
be consistent with the predictions from the corresponding mean-field rate equations. If spreading occurs 
only through nearest-neighbor hopping, small population clusters emerge; yet the requirement of an in- 
termediate species cluster obviously disrupts spatio-temporal correlations between predator and prey, and 
correspondingly eliminates many of the intriguing fluctuation phenomena that characterize the stochastic 
spatial LV system. 

PACS. XX.XX.XX No PACS code given 

J> ■ 1 Introduction [TU], and have been invoked to capture the remarkable 

coexistence of three cyclically competing types of Califor- 
. Spatially extended evolutionary game theory models have nian lizards J23.J.241, and of three strains of E. coli bacteria 
J^' found widespread applications in population dynamics, in microbial experiments gS]. It was demonstrated that 
■ - ■ with the aim to explore the complex spatio-temporal pat- stochastic lattice versions of both LV and RPS models 
terns that characterize ecological and biological systems may exhibit spatially correlated evolving patterns which 
in which multiple species interact cooperatively or com- turn out to be essential in maintaining species coexistence 
petitively [1 , 2 , 3 , 4 , 5 .EUZIISIE] • Various simplified models over extended time periods [ni[Igl[Il[Igl[Tl[gTl[ggl[gg[g71 
have been employed to qualitatively capture and quan- [55]. Ecological systems in nature of course typically con- 
titatively understand the primary features of interacting tain many more interacting species; it is therefore vitally 
multi-species systems. Two prominent non-trivial paradig- important to understand if the intriguing observed fea- 
matic representatives are the two-species Lotka-Volterra tures in two- and three-species models are generic to multi- 
model (LV) [Tni[IIl[Tl|Tl[Tl[TSl[Tni[l71ll81[Il] that describes species systems as well, and under which circumstances, 
predator-prey competition and coexistence, and the three- if at all, an effective description in terms of fewer species 
species rock-paper-scissors (RPS) model [Il[l|3,4,20,21, is appropriate. In this paper, we consequently aim to ex- 
that encapsulates cyclic competition through preda- plore the relationship between both cyclic and hierarchical 
tion. LV and RPS model variants effectively mimic the three-species predator prey systems and the two-species 
dynamical evolution of coexisting species in interacting LV model. To this end, we shall consider both well-mixed 
multi-species systems or autocatalytic chemical reactions 
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settings that are properly described by the correspond- 
ing mean-field rate equations, and stochastic dynamics on 
a two-dimensional square lattice with periodic boundary 
conditions, where population spreading occurs through 
nearest-neighbor hopping and/or particle exchange. 

The classic LV model captures the dynamical evolu- 
tion of two competing species of predators and prey; in the 
well-mixed mean- field limit, the ensuing coupled determi- 
nistic nonlinear ordinary differential equations yield regu- 
lar periodic population oscillations. However, this original 
non-spatial LV model is often criticized for its oversim- 
plification [5] and lack of robustness: it is mathematically 
unstable against various model modifications and varia- 
tions [29]. Therefore, hopefully more realistic stochastic 
spatial LV model extensions have been studied extensively 
both analytically and computationally; this includes the 
stochastic lattice-based LV model with unrestricted site 
occupancies [151150] and model versions with restricted lo- 
cal carrying capacity, where only a limited number of indi- 
viduals may occupy each lattice site [T2l[T5l[T6lll9ll31ll32] . 
Systems with restricted site occupancy have been shown 
to display a critical extinction threshold for the predator 
population, where in the thermodynamic limit a contin- 
uous phase transition takes place from an active coexis- 
tence state to an inactive absorbing state, governed by the 
power laws of critical directed percolation [T^fT^ITB II 1 9 lIMl 
[521I5?1[M1I^ ■ In the active (quasi-)stationary coexistence 
state, stochastic spatial LV models display essential and 
robust features such as species clustering and the emer- 
gence of correlated spreading activity fronts [T^ITSllTSlfroi 
[36ll37ll38ll39ll40] . which are associated with the presence 
of persistent erratic population oscillations [T^IT^ITSlfTBl 
[T5iri9,32.35.41 : in spatial stochastic systems, these char- 
acteristics even pertain when the original LV predation 
reaction is split up into two independent processes [42^ . In 
addition, fitness enhancement for both species is induced 
by significant spatial variability in the predation rates that 
control the interactions between predators and prey [30,. 

Cyclic competition in three-species RPS systems has 
been suggested to provide a robust mechanism to pro- 
mote species coexistence. For non-spatial stochastic RPS 
systems, the temporal evolution always reaches one of 
three extinction states in which two of three species disap- 
pear [inmilllS] ■ In contrast, sufficiently large spatially ex- 
tended RPS models, particularly stochastic two-dimensional 
RPS systems, are characterized by the coexistence of all 
three particle species. When the cyclic competition reac- 
tions do not conserve the total population density, com- 
plex spatio-temporal structures such as spiral patterns 
emerge in two-dimensional RPS systems [2^l?7ll^ll46IIT7l 
H51l49 , 50 . Yet in the simplest spatial stochastic RPS real- 
izations of cyclic three-species competition, a conservation 
law for the total number of reactants is incorporated, and 
no spiral patterns become apparent in lattice simulations, 
but instead fluctuating species clusters are observed [CTl 
[521155] . However, even though spatially correlated clusters 
emerge in RPS systems with conservation law, quenched 
spatial disorder in the reaction rates does not evidently in- 
fluence the dynamical evolution, provided these rates for 



the different species remain comparable [52] , which locates 
the coexistence state far away from the "corners" of con- 
figuration space. In this present work, we shall therefore 
study the effect of quenched spatial disorder on RPS sys- 
tems with strongly asymmetric reaction rates, whose co- 
existence states are shifted closely towards the "corners" 
of configuration space. 

Since both the two-species LV model and cyclic three- 
species RPS systems have attracted considerable attention 
in the literature, understanding the emergence of stable 
biodiversity in population dynamics also requires a better 
grasp on potential connections between these two paradig- 
matic model systems, which might in turn perhaps enable 
us to reduce complex interacting multi-species systems to 
effective models with fewer degrees of freedom. Now let 
N represent the typical size of the "living space" for each 
species (e.g., in a lattice-based simulation, N is the typi- 
cal size of the lattice) , and consider a cyclic three-species 
RPS model with interaction rates chosen in such a manner 
that one species reaches a persistently large population 
density, say of 0(1), in the (quasi-)steady state (hence- 
forth referred to as the majority species), while the other 
two population densities are reduced to 0{1/N) (minor- 
ity species). Through treating the majority species in the 
system in effect equivalent to the empty states in a two- 
species LV model, the predation processes between the 
majority and the two minority species in such a "corner" 
RPS system can be rewritten as the elimination and repro- 
duction processes for the predators and prey, respectively, 
in the LV model. That is, it should be possible to ap- 
proximate the evolutionary dynamics of the two minority 
species in the "corner" RPS model by the LV model; and 
the error limit of this approximation has been shown to be 
0(l/iV^) [53]. Another realization of an interacting three- 
species system is a hierarchical "food chain" model, which 
is generated by inserting one intermediate species between 
the predator and prey in two-species LV model; the preda- 
tions between predators and prey are then only indirectly 
generated via their interactions with the intermediate pop- 
ulation. It is then an interesting question whether one may 
still employ the two-species LV model with fewer degrees 
of freedom, to effectively approximate the full stochastic 
dynamic behavior of the predator and prey populations 
in such spatially extended three-species "food chain" sys- 
tems. 

The goal of this work is to explore these two distinct 
ideas as to how a stochastic spatial predator-prey system 
with originally three interacting species might be mapped 
onto an effective two-species LV model, and thereby bet- 
ter comprehend potential mechanisms for complexity re- 
duction in multi-species population dynamics. Our main 
results can be summarized as follows: 

(1) We demonstrate that the two minority species in 
three-species "corner" RPS system exhibit approximately 
the same evolutionary dynamics as the predators and prey 
in the two-species LV model: We observe quantitatively 
similar population oscillations and asymptotic densities 
in the (quasi-)steady state, noticeable species clustering, 
and robust fitness enhancement due to spatial variabil- 
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ity in the predation rate, i.e., larger population densities 
for both species, associated with reduced relaxation times 
towards coexistence, and more localized species clusters. 
These features for the two minority species in the "cor- 
ner" RPS model closely resemble previous findings for 
the efi^ect of quenched spatial disorder in the two-species 
LV system '30]. We conclude that these numerical results 
confirm the assertion, based on the mean-field analysis, 
that two-species LV models satisfactorily approximate the 
stochastic spatial evolutionary dynamics of the two minor- 
ity species in the three-species "corner" RPS system. 

(2) We also study the hierarchical three-species "food 
chain" model. In its spatially well-mixed version, there ex- 
ists an extinction threshold for the predator species, just 
as in the two-species LV model with restricted site oc- 
cupation. We find that nearest-neighbor particle pair ex- 
change processes alone are strong enough to wash out the 
emergence of species clusters, eliminate the influence of 
quenched spatial disorder and thus promote the appear- 
ance of mean-field like behavior. Our numerical simula- 
tion results for the coexistence states of these effectively 
well-mixed systems are consistent with the mean-field pre- 
dictions. In contrast, if species spreading happens exclu- 
sively through nearest-neighbor hopping, the spatial "food 
chain" system does not behave according to the predic- 
tions from the mean-field approximation, and rather dis- 
tinctive species clusters emerge as observed in two-species 
LV systems. However, due to the necessity of interspersed 
clusters of the intermediate species, the effective reac- 
tion boundaries between the predator and prey clusters 
in these spatially segregated systems are strongly sup- 
pressed, and correspondingly spatial disorder in the re- 
action rates still displays only a minute influence on the 
evolutionary dynamics of "food chain" model. That is, the 
predator and prey in the hierarchical three-species "food 
chain" model in a stochastic spatial setting do not behave 
akin to their LV counterparts. 

The organization of this paper is as follows: In Sec. [U 
we deflne the lattice-based stochastic rock-paper-scissors 
(RPS) model with conserved total particle number, as 
well as the stochastic Lotka-Volterra (LV) model with 
site occupation number restriction, briefly discuss well- 
established results from mean-field theory, and analyze 
the quantitative relationship between the (quasi-)steady 
states of the "corner" RPS model with strongly asym- 
metric rates and its associated two-species LV model. In 
Sec. [21 we explain how we implement our Monte Carlo 
simulations on a two-dimensional lattice with restricted 
site occupancy, allowing at most a single individual of ei- 
ther species per site, and introduce the relevant physical 
quantities of interest. Then we present and analyze our 
simulation results for the corner RPS and associated LV 
models, both in the absence and presence of quenched re- 
action rate disorder, in Sec. [31 In Sec. [H we describe the 
hierarchical three-species "food chain" model and discuss 
the corresponding analytic results from the mean-field ap- 
proach. Then the numerical results based on Monte Carlo 
simulations on a two-dimensional lattice with at most one 



individual per site are presented. Finally, Sec. [SI concludes 
with a discussion and interpretation of our findings. 

2 Strongly asymmetric "corner" RPS model 
and mean-field analysis 

The rock-paper-scissors model (RPS) consists of three 
zero-sum predator-prey interactions [2]. Let A, B, and C 
represent the three interacting particle species; the RPS 
model is then described by the following binary reactions: 

A + B ^ A + A with rate A ; 

B + C ^ B + B with rate a ; (1) 

C + A ^ C + C with rate fj, 

(here, we use A, a, to label the corresponding reaction 
rates to resemble the notation used in the two-species 
LV model [T51ITM5D] ). It is important to note that the 
total particle population number is always conserved by 
the reactions ([T]). In addition, in our lattice Monte Carlo 
simulations we shall allow for nearest-neighbor exchange 
processes X + Y ^ Y + X with X,Y e {A, B, C). We 
also impose the constraint that at most a single particle 
(of either species) is allowed on each lattice site to mimic 
limited local carrying capacities that result from limited 
natural resources. In this present study, the total popu- 
lation density is always set to 1, and the lattice is hence 
fully occupied. 

Let a(t), b{t), and c{t) represent the (spatially av- 
eraged, see Eq. ^ below) population densities or con- 
centrations of species A, B, and C, respectively. Within 
the mean-field approximation, the associated coupled rate 
equations 

dta{t) ^ait)[Xb{t)~fic{t)] , 
dtKt)^b{t)[acit)-Xa{t)] , (2) 
dtc{t) = c{t) [na{t) -<Tb{t)] 

yield one marginally stable reactive fixed point, describing 
a three-species coexistence state (a*,&*,c*) = p{a,fj.,X)/ 
{a + fi + X), where p = a{t) -\- b{t) -\- c{t) denotes the 
conserved total particle density [5T[[52] (with p = 1 in 
this study). The three absorbing states (p, 0,0), (0, p, 0), 
and (0, 0, p) are all unstable in the mean-field approxima- 
tion, but one of them will be eventually reached in any 
finite stochastic system, after a characteristic extinction 
time that increases exponentially with system size [26ll28l 
152] . In our simulations, we employ sufficiently large lat- 
tices that extinction events are extremely unlikely within 
the run times. Therefore, systems with comparable reac- 
tion rates X ^ a ^ ^ approach a (quasi-)steady state 
far away from the "corners" of configuration space. It has 
been demonstrated that stochastic fluctuations are com- 
paratively small in such systems on two-dimensional lat- 
tices, and quenched spatial disorder in the rates has only 
minor influence on their dynamical evolution ,52 . 

However, in the strongly asymmetric limit A 3> cr, ^, 
the mean-fleld stationary densities become to leading or- 
der in 1/A: (a*,6*,c*) « /o(f , ^,1 - and the corre- 
sponding flxed point moves to the vicinity of one of the 
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corners of configuration space. Due to the relatively large 
reaction rate A, species C acquires a very high popula- 
tion density (almost saturated, i.e., c* « p) in the (quasi- 
)steady "corner" state. Since the C particles almost uni- 
formly fill the lattice, the minority species A / B will es- 
sentially always encounter a nearest-neighbor partner C to 
undergo the third / second reaction in ([1]). The evolution 
of the rare species A and B in the RPS corner state can 
therefore be approximated by the following two-species 
Lotka-Volterra model reactions: 



Its eigenvalues are 



A — > with rate p. = pc* 

A + B ^ A + A with rate A ; 

B ^ B + B with rate a = a c* 



(3) 



Here, the symbol denotes the empty state. Notice that in 
this stochastic effective two-species LV model on a lattice, 
empty sites inevitably appear in the system, no matter 
whether the lattice is initially fully occupied or not. In 
effect, the original exchange processes of minority particles 
A 01 B with the majority C species are replaced with 
nearest-neighbor hopping to empty sites in the resulting 
two-species LV model. 

In the mean- field approximation, the LV model (jSj, 
with the constraint that at most one particle is permitted 
to reside on each site, is described by the following two 
coupled rate equations, where the total density p assumes 
the role of the overall carrying capacity [roil29) : 



dta{t) = a{t) [Xb{t) - p] 



Xa{t)b{t) . (4) 



These also follow directly from the RPS rate equations 
([2]) by setting c{t) « c*; more precisely, in the second 
rate equation for b(t), employing the exact conservation 
relation c{t) = p — a{t) — b{t) properly incorporates the 
limited local carrying capacities, and a ^ a p. The sta- 
tionary states of the coupled rate equations @ consist 
of the linearly unstable extinction state (a, 6) = (0,0), 
another absorbing state (0,p) that is linearly stable for 
small X < Xc — p/p ~ p, and the coexistence state 

{as,bs) — (^ ^^xp i^P^ f) ' ^) ^^^^ exists and becomes 
linearly stable if A > Ac- Therefore, in the thermodynamic 
limit the system displays a continuous non-equilibrium 
phase transition from an inactive, absorbing state to the 
active coexistence state at the critical predation rate Ac. 
The universal power laws that emerge near the extinction 
threshold are characterized by the critical exponents of 
directed percolation P^[TC1[T^[55] . 

Linearizing around the active coexistence fixed point 
leads to 



dt 6a 
dtSb 



(5) 



where Sa{t) ~ a{t) — and 6b(t) ~ b{t) — b^, and with 
the linear stability matrix 



L = 



1 







Xpa-{Xp — p) 



Xp{Xp + &) \-p{Xp + ay -ap{Xp + a) 



(6) 



e± = -aA(2Ap)-i[l ± ^1 - 4Apa-i(A/9/i-i - 1)] ■ (7) 

In the limit X ^ a, p in the two-species LV model ([3]), 
we have Xpp^^ = Xp{pc*)~^ ~ Xp~^ 3> 1 and similarly, 
Xpa~^ S> 1, whence the eigenvalues e± for the active coex- 
istence fixed point turn into a complex conjugate pair with 
negative real part. This demonstrates that, in the strongly 
asymmetric limit X ^ a, p, the nature of the nontrivial 
coexistence fixed point is that of a linearly stable spiral 
singularity, and the (quasi-)steady state consequently is 
approached in an exponentially damped oscillatory man- 
ner. From the imaginary part of the complex conjugate 
pair ([7]), we infer the characteristic oscillation frequency 
/ = w/(27r)« V^/(27r). 

Moreover, there are two points worthwhile noticing: 

(i) The strong asymmetry condition A ^ p,a is the pre- 
requisite for validating the approximative representation 
of the stochastic model ^ through the effective model ([3]); 
therefore, in the resulting two-species LV system Q the 
predation rate A is always larger than the critical threshold 
Ac — pc* / p ^ A, and correspondingly the system resides 
deep in the active coexistence state. 

(ii) When implementing the Monte Carlo algorithm for 
the effective model ([3]), instead of only introducing a lo- 
cal growth limit for the prey species B (which in the 
mean-field limit would be described by the rate equation 
dtb{t) = ab{t)[l - b{t)/p] - Xa{t)b{t) ^MSM): we im- 
pose a maximum total population carrying capacity on 
each site in the simulation. Deep in the coexistence state, 
one would expect either variant to induce similar evolu- 
tionary dynamics for the model First, since A 3> cr, 
the stationary mean-field coexistence state concentrations 
for both model variants, namely with growth-limiting con- 
straint on either the total or just the prey population, 

are essentially identical, i.e., {as,bs) ~ , 

Second, we shall see that in our stochastic spatial sim- 
ulations, the characteristic oscillation frequencies, which 
are inferred from the imaginary part of the eigenvalues 
e± in the nontrivial deep coexistence state, remain basi- 
cally unchanged even if we only impose spatial occupancy 
restriction on the prey population [TSlfTO] . 

Altogether, it thus appears reasonable and legitimate 
on the ground of both mean-field theory and heuristic 
considerations that the stochastic spatial "corner" three- 
species RPS system approaches the two-species Lotka- 
Volterra model. In the following section, we shall demon- 
strate through detailed Monte Carlo simulations that this 
claim is valid. 



3 Monte Carlo simulation results for the 
corner RPS and the associated LV models 

We implement stochastic individual-based Monte Carlo 
simulations for both the "corner" RPS and the associated 
LV model on two-dimensional square lattices (typically 
with N = 256 x 256 sites) with periodic boundary condi- 
tions according to the reaction schemes ((Ij and ^ that 
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define these systems. The RPS model is moreover subject 
to nearest- neighbor exchange processes X + Y ^ Y + X, 
where X,Y G {A,B,C), while in the LV system we also 
allow nearest-neighbor hopping X + ^ ^ ^ + X, in addi- 
tion to particle exchange. In either situation, a maximum 
occupancy of a single particle (of either species) is im- 
posed on each lattice site to mimic a finite local carrying 
capacity resulting from limited natural resources. At each 
time step, one individual of any species is chosen at ran- 
dom, and simultaneously one of its four nearest neighbors 
(which might be an empty site in the LV system) is se- 
lected randomly. Subsequently, the center individual and 
its chosen neighbor undergo the various reactions defined 
in the models according to the respective associated rates; 
otherwise, exchange (or hopping for LV) processes take 
place between the center particle and its chosen neighbor. 
Once on average all individual particles on the lattice have 
had a chance to be selected as center individual, one Monte 
Carlo step (mcs) is completed; therefore, the infinitesimal 
simulation time step is 6t ~ M~^, where M is the to- 
tal number of individuals on the lattice present at that 
time (a fixed constant only for the RPS system). In the 
following, we are going to investigate four model variants: 

1. "RPS" — "Corner" rock-paper-scissors model with 
strongly asymmetric rates, A ^ tr, /i, whence the system 
is thus going to eventually reach one of the corners of con- 
figuration space. In addition, we take the lattice for this 
model to always be fully occupied, p — p{0) — 1; therefore, 
only exchange processes take place in the simulation. 

2. "RPS with disorder" — In order to study the effect 
of quenched spatial disorder, we will treat the predation 
rate A at each lattice site as a random number drawn from 
a normalized Gaussian distribution truncated at one stan- 
dard deviation on both sides: For example, A ^ M{m, n) 
implies that the predation rate A is drawn from a trun- 
cated normal distribution on the interval [m — n, to -1- n], 
centered at the value to with standard deviation n < m. 
In the simulation, the randomized rates are attached to 
each site (rather than to each individual), and held fixed 
(quenched) for a given disorder realization: The predation 
rates A on all sites are determined at the beginning of 
each single Monte Carlo run, and remain unchanged until 
the next single run is initiated. Again, the total popula- 
tion density in the simulation is always set to p — p(0) = 1 
and thus only exchange processes may take place. All other 
rates except A stay the same as in the "RPS" models. 

3. "LV" — Lotka-Volterra model with rates that orig- 
inate from the rates in the corresponding corner RPS sys- 
tem: If rates A, a, p are used in the RPS model, in order to 
appropriately compare the numerical results for both sys- 
tems, we take A,(T = ac*, fi — pLc* for the associated LV 
model, where c* represents the mean-field population den- 
sity for the majority species C in the (quasi-)steady state 
in the asymmetric corner RPS model, as determined from 
the stationary solutions to the rate equations Due to 
the absence of a conservation law for the total population 
density, empty sites inevitably occur in the course of the 
simulation runs, whence nearest-neighbor hopping is al- 



lowed in the simulations in addition to particle exchange 
processes. 

4. "LV with disorder" — We treat the predation rate A 
at each lattice as a quenched random variable drawn from 
the same truncated normalized Gaussian distribution as 
in "RPS with disorder" . Meanwhile, all other rates remain 
the same as those in the "LV" models, and again both 
exchange and hopping processes may take place. 

Typically, in order to study the connection between 
the corner RPS and its associated effective LV model, 
we shall numerically investigate their emerging spatio- 
temporal structures through instantaneous snapshots of 
the particle distributions on the lattice, the temporal evo- 
lution of the (spatially averaged) population densities, e.g. 
for species A: 

a{t)^{nA{j,t))^^Y.'^A{j,t) , (8) 

the temporal Fourier transform of these signals, 

«(/) = j ait) e^-/* dt , (9) 

and the equal-time two-point correlation functions in the 
(quasi-)steady state, for example 

CAB{x,t) ^ {haU + x,t)nBij,t)) - a{t)b{t) , (10) 

where j and j + x represent square lattice site indices. 

3.1 Self-organization in the coexistence phase of the 
two minority species in the corner RPS model 

In Fig. Hal the population densities for the minority species 
A and B in an asymmetric "corner" RPS system model on 
a square lattice with rates A = 0.8 and a = p = 0.07 are 
plotted as functions of simulation time t (in mcs). While 
the population densities for species A and B as obtained 
from the mean-field approximation for the RPS model 
would be identical, j+f:p^(cr, p) w (0.074, 0.074), the 
corresponding population densities in the associated LV 
model assume unequal values, {as,bs) = (0.065,0.075) in 
the (quasi-)steady state. The LV predictions are remark- 
ably close to the measured values (0.064 ± 0.003, 0.074 ± 
0.002) (averaged over 50 runs a.tt = 1000 mcs, see Fig. [Ta|) . 
clearly distinct for species A and B. In particular, the per- 
sistent unequal population densities (starting from t ~ 600 
mcs in Fig. [Ta|) support the population density inequality 
for minority species A and B in the (quasi-)steady state as 
statistically reliable. At the beginning of the simulations, 
we observe marked population oscillations, see Fig.[Tal re- 
flecting the formation of complex spatio-temporal struc- 
tures, i.e., expanding activity fronts as familiar in the 
stochastic spatial LV system, noticeable as minority speci- 
es clusters in Fig. Ilbl The emergence of such localized 
species clusters promotes species coexistence, reduces the 
oscillation amplitudes in the system, and eventually allows 
the system to evolve into the (quasi-)steady coexistence 
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Fig. 1: (Color online.) (a) Temporal evolution for the population densities of the minority species A and B for an 
asymmetric RPS system with rates A — 0.8 and a = fj, = 0.07 on a square lattice with N = 256 x 256 sites, periodic 
boundary conditions, and initial densities a(0) = 6(0) = 0.05, c(0) — 0.9; the data are averaged over 50 simulation 
runs, (b) Corresponding snapshot of the spatial particle distribution aX t = 1000 mcs (yellow/light gray: A, red/gray: 
B, black: C); note that there is no empty site in the lattice, since p = a(0) -|- 6(0) -|- c(0) = 1 is conserved under the 
temporal evolution. 



state. In light of these observations, we conclude that the 
temporal evolution of the two minority species A and B 
in the stochastic spatial corner RPS model with strongly 
asymmetric rates can indeed be faithfully approximated 
by the predator-prey behavior in the related two-species 
LV model. We shall further support this assertion through 
additional observations in the following subsection. 



3.2 Fitness enhancement of the minority species due 
to spatial variability in the corner RPS model 

In order to further and quantitatively characterize the con- 
nection between the strongly asymmetric corner RPS sys- 
tem and its corresponding effective LV model, we next 
study the influence of quenched spatial disorder on the 
evolution of either system. Since spatial variability turns 
out to have very similar effects on species A and B, we 
mainly show the results associated with species A; a com- 
plete listing of data for characteristic observables is pro- 
vided in Table [T] We depict the temporal evolution of the 
population density a{t) = {nA{j,t)), its Fourier transform 
a(/), and the spatial auto- and cross-correlation functions 
Caa and Cab, respectively, in the (quasi-)steady state 
(here, obtained at t = 3000 mcs) in Fig. [21 The default 
rates are set to A = 0.5, a = fi — 0.03 in the corner RPS 
model, and correspondingly A = 0.5, a — p, — 0.027 in the 
associated LV model, and A ~ J\f {0.5, 0.4) in both corner 
RPS and LV systems with quenched spatial disorder in 
the rate A. Note that within the mean-field approxima- 
tion, the population density for species A in the (quasi-) 
steady state is as ~ 0.048, and the characteristic critical 
predation rate threshold is Ac = 0.027 < A, guarantee- 



ing that the associated LV system is located in the active 
(quasi-)steady coexistence region of the phase diagram; it 
is also worth noticing that the minimum predation rate in 
the model variants with quenched spatial disorder is 0.1, 
which is still well above the critical threshold, eliminat- 
ing the possibility of any species extinction in our large 
lattices, even in the presence of spatial rate variability. 

Figure [2a| shows the simulated average population den- 
sities as function of time for the minority species A in 
both corner RPS and associated LV systems, in the ab- 
sence and presence of quenched spatial disorder. In the 
stationary state, the A species concentrations (averaged 
over 50 runs at t = 3000 mcs) become 0.057 for the corner 
RPS model, 0.060 for the associated LV model, 0.063 for 
the RPS model with quenched spatial disorder, and 0.068 
for LV model with spatial disorder (see Table [T]). Thus, 
within the error bars, the population density (0.057) for 
species A in the "pure" corner RPS model coincides with 
the simulation result (0.060) for the associated stochastic 
spatial LV model. Moreover, we observe that spatial vari- 
ability in the corner RPS model enhances the fitness of the 
minority species, as measured by the asymptotic station- 
ary concentrations, by ~ 11%. This fitness enhancement 
is in qualitative accord with earlier investigations on the 
effect of quenched randomness in the predation rate for 
the the two-species LV model without site restriction [5D] , 
but turns out to be quantitatively less significant than 
the ^ 25% found there; it is also slightly smaller than the 
^ 13% enhancement obtained directly in our simulations 
of the two-species LV model with restricted site occupa- 
tion number, see Fig. [2al and Table [1] 

We may attribute these differences to the fact that spa- 
tial disorder also affects the stationary population density 
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(a) Temporal evolution a{t). 
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(c) Correlation function Caa{x) 



(d) Cross-correlation Cab{x) 



Fig. 2: (Color online.) Quantitative observables for a stochastic comer RPS model with rates A — 0.5, a = fi — 0.03, 
and initial densities a(0) = 6(0) = 0.02, c(0) = 0.96, and its corresponding LV system with A = 0.5, a = jj. = 0.027 and 
a(0) = 6(0) = 0.03 from simulation runs on a square lattice with N = 256 x 256 sites (periodic boundary conditions). In 
the presence of spatial disorder, A ^ Af{Q.5, 0.4). All simulation results were averaged over 50 runs, and the correlation 
functions in (c) and (d) were measured at t — 3000 mcs. 



c* of the majority species C in the (quasi-)steady state of 
the corner RPS model, which in turn implies locally vary- 
ing effective rates a and fi for the associated two-species 
LV model. This becomes apparent already on the mean- 
field description level, where c* — x+a+n ^^'^ ^ ~ ^ 
jl = /j,c* . Then for the minority species concentration deep 
in the coexistence state of the corner RPS model one ob- 
tains Os ~ J = x+f+M' Simil^^ly: for the two-species LV 
model with site restrictions, as ~ XT^' "whereas = j 
in the LV model with infinite local carrying capacity. In 
the latter broad distribution of predation rates 

A strongly biases the system towards small values that 
considerably enhance the associated (quasi-)steady state 
density of A species [22 . For the RPS system or the LV 
model subject to site occupation restrictions, the addi- 
tional rates in the denominators reduce the significance of 
this low-A bias and ensuing concentration enhancement. 



and consequently the overall fitness enhancement is less 
pronounced. 

As one would expect, the emergence of species clus- 
ters results in remarkable population oscillations at the 
beginning of the simulation runs, as seen in Fig. I2al How- 
ever, as shown in Fig. I2al the transient period for the 
system evolving from the initial to the eventual active 
(quasi-)steady state is quite distinct for the different mod- 
els. The corresponding Fourier-transformed density signal, 
which reflects both the characteristic oscillation frequency 
in the steady state and the associated characteristic relax- 
ation time, is displayed in Fig. I2bl The Monte Carlo sim- 
ulations exhibit the same peak frequency in the (quasi-) 
steady state for all four models, / ~ 0.0014; this charac- 
teristic peak frequency in the stochastic lattice systems is 
markedly reduced by ~ 60% as compared with the (lin- 
earized) mean-field prediction / « \//Ia-/(27r) « 0.0043. 
This remarkable decrease in the typical population oscil- 
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as 


bs 


ta/b 


Iaa 


Ibb 


Iab 


RPS model 


0.057 ± 0.003 


0.060 ± 0.002 


~ 1800 mcs 


13.2 ± 0.6 


16.0 ± 2.0 


~ 13.0 


RPS with disorder 


0.063 ± 0.002 


0.067 ± 0.003 


~ 800 mcs 


9.4 ± 0.2 


13.0 ± 1.0 


~ 8.0 


LV model 


0.060 ± 0.003 


0.063 ± 0.003 


~ 2000 mcs 


12.8 ± 0.5 


15.5 ± 0.5 


~ 13.0 


LV with disorder 


0.068 ± 0.002 


0.070 ± 0.003 


~ 1300 mcs 


10.3 ± 0.3 


11.2 ± 0.2 


~ 9.0 



Table 1: Characteristic physical quantities measured by Monte Carlo simulations for the four stochastic spatial RPS 
and LV models whose features are displayed in Fig. [5] The population densities {Us, bg) in the (quasi-)steady state are 
obtained as the averages over 50 simulation runs a.t t = 3000 mcs (see Fig. I^al) . the associated errors are measured 
as the standard deviations of the 50 data; and the relaxation time r is found as the inverse of the full width at half 
maximum of the frequency peak in Fig. I2bl The correlation lengths Iaa and Is b are extracted by fitting the simulation 
data for Caa and Cbb to exponential functions (see Fig. I2cp : and the separation length Iab is measured as the site 
location with the maximum value for Cab (see Fig. I2d[) . 



lation frequency can be attributed to the renormalization 
by stochastic fluctuations in the spatial system 18, 19,55 . 
Moreover, the relaxation time r = 1/^, where F denotes 
the full width at half maximum of the oscillation peak in 
Fourier space, characterizes the typical time scale for the 
system to relax towards the (quasi-)steady state, and thus 
represents a measure of stability for these systems against 
external perturbations. Significantly, as demonstrated by 
the data in Fig. [2b] and Table [1] the relaxation time in the 
RPS model with quenched spatial disorder is reduced by 
^ 50% to ^ 800 mcs as compared with t ~ 1800 mcs in 
the corner RPS model without spatial rate variability. 

The enhancement of both predator and prey fitness re- 
sulting from spatial variability in the stochastic lattice LV 
model is based on the emergence of more localized cor- 
related spatio-temporal structures [50]. In all four model 
systems investigated here, such species clustering appears 
as well, see Fig. [Tb| Compared to the RPS/LV model 
without spatial disorder, the corresponding RPS/LV sys- 
tem with spatial variability in the reaction rate A displays 
larger initial oscillation amplitudes (c.f. Fig. Eaj) . imply- 
ing that spatial disorder in the predation rate tends to 
gather species closer, resulting in more localized and thus 
enhanced population growth spurts. To better understand 
the underlying spatial structures in our systems, we depict 
the static two-point auto-correlation functions in Fig. \Tc\ 
and cross-correlation functions in the (quasi-) steady state 
in Fig. I2dl From these plots, we extract the typical cor- 
relation length Iaa-, which measures the spatial extent of 
species A clusters by fitting the simulation data Caa with 
an exponential function exp(— and obtain the 
predator-prey separation length Iab-, defined as the lo- 
cation where Cab assumes its maximum positive value. 
In the corner RPS model, the correlation and separation 
lengths turn out to be Iaa = 13.2 (in units of lattice spac- 
ing) and Iab = 13.0, respectively, both of which coincide 
within the statistical error bars with the simulation results 
for the associated LV model, for which we find Iaa — 12.8 
and Iab = 13.0, as listed in Table [T] However, in the cor- 
responding disordered RPS system, we measure markedly 
reduced values for the correlation and separation lengths, 
namely 9.4 and 8.0. This closer clustering near sites with 
locally favorable rates permits an overall larger number 
of population centers, and consequently enhances the net 



population densities for the minority species in the system 
with spatial variability in the rate A. 

In summary, as suggested on the basis of mean-field 
considerations, we find that the minority species A and 
B in our spatially extended stochastic RPS model with 
strongly asymmetric rates behave like the predator and 
prey in the coexistence phase of the two-species Lotka- 
Volterra model. Indeed, similar to our earlier findings for 
the lattice LV model [30], we observe that quenched spa- 
tial disorder can markedly enhance the fitness of both mi- 
nority species. This feature distinguishes the corner RPS 
model from the perhaps more common RPS system with 
comparable reaction rates, whose stationary state resides 
close to the center of configuration space; in that situation, 
quenched spatial reaction rate disorder hardly influences 
the dynamical evolution of either species [521150] . 



4 Hierarchical three-species "food chain" 
model 

In Section [3j we have demonstrated how the properties 
of the coexistence state in the two-species Lotka-Volterra 
model emerge as the rock-paper-scissors system is moved 
towards one of the corners of configuration space by choos- 
ing strongly asymmetric rates. Another natural idea to 
generate the LV model is to consider a hierarchical three- 
species "food chain" system, where an intermediate species 
C is inserted between the predators A and prey B. The 
question then becomes: will a stochastic hierarchical three- 
species food chain model on a lattice in a certain limit 
again recover the properties of the spatial two-species LV 
system? To this end, we investigate the following cou- 
pled stochastic reactions that define our three-species food 
chain: 

A — ^ with rate fi 

A + C ^ A + A with rate A 
C + B ^ C + C with rate A; (11) 

B ^ B + B with rate a . 

Species A and B in this model behave as predators and 
prey, respectively, while the intermediate population C 
preys on species B and is preyed upon by species A. In 
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order to perhaps closely approach the two-species LV sys- 
tem and mainly study the dynamic behavior of species A 
and B, we here use identical reaction rates A for the A-C 
and B-C predation processes. We remark that upon iden- 
tifying the empty state with a fourth species, this system 
can be mapped onto a cyclic four-species predator-prey 
model, see, e.g., Refs. t56''57', with a peculiar choice of 
rates. In our Monte Carlo simulations for a spatially ex- 
tended system on a square lattice with periodic boundary 
conditions, we allow at most one particle per site. The 
resulting mean-field rate equations now read: 



dta{t) = a{t) [Xc{t)-^l] , 
dtc{t) = Xcit)[b{t)-a{t)] , 
dtb{t) = ab{t) [1 - a{t) - b{t) 



(12) 



c{t)] - Xb{t)c{t) 



since only the reproduction process for the species B re- 
quires the presence of an empty site in its immediate 
neighborhood, whereas newly generated A and C parti- 
cles just replace their respective prey. As stationary so- 
lutions of the coupled rate equations P^ . one finds one 
linearly unstable absorbing point (a, c, fe) — (0,0,1), an 
entire absorbing line (0,c, 0), where c G [0,1], and one 
active coexistence fixed point: 



{as,Cs,bs) = 



/Act — /i(A + (j) n Act — /i(A + ct) 



2Act 



' A 



2Act 



.(13) 



As in the two-species LV model with site occupation re- 
strictions, there exists a critical threshold for the preda- 
tion rate Ac = in this hierarchical food chain model, 
where in the thermodynamic limit a phase transition oc- 
curs from the inactive absorbing states to the active coex- 
istence state: When the predation rate is below the thresh- 
old (A < Ac), the system evolves towards extinction for the 
A and B population (0, c, 0); otherwise, the system eventu- 
ally reaches the three-species coexistence state (a^, Cs,bs). 
However, in stark contrast with the mean-field analysis 
for the two-species LV model with site restrictions, the 
population densities for the predators A and the prey B 
in this food chain system are always the same in the sta- 
tionary coexistence state, in order for the intermediate 
species C to possess a nonzero stationary population den- 
sity. Other important distinct spatio-temporal properties 
for the spatially extended hierarchical food chain model 
will be discussed in the following two subsections. 



4.1 Monte Carlo simulations: mean-field like behavior 

In the Monte Carlo simulation runs, we first permit both 
nearest-neighbor particle exchange and hopping processes 
to generate population spreading. In Fig. l3a[ we depict the 
typical normalized trajectories for the population densities 
in systems with varying predation rates A. The other rates 
jjL = 0.01, CT = 0.1 are unchanged, whence the critical pre- 
dation rate threshold value remains fixed at Ac = ~ 
0.011. When the predation rate is below this threshold, the 
trajectories (A = 0.0001, 0.005 in Fig. [5a]) reach the bound- 
aries in the phase space spanned by the the population 



densities, which constitute inactive absorbing states. Upon 
increasing the predation rate (A = 0.02,0.1 in Fig. [3a|) . 
we observe that the system evolves to the active coexis- 
tence state. While for small predation rates just above the 
threshold (A = 0.02 in Fig. [5a|) the fixed point is a node, 
it changes character and becomes a spiral focus for larger 
values, as evidenced by the spiral trajectory in the phase 
graph (A = 0.1 in Fig. [5a|) . This fixed point evolution 
as function of the predation rate is remarkably similar to 
the situation in two-species LV model with site restric- 
tion [TS]. As discussed earlier, in the active coexistence 
state the predator A and prey B species should possess 
the same population densities in order to let the interme- 
diate species C also have nonzero stationary population 
density (see the dashed line in Fig. [5a|) . However, when A 
becomes large relative to the rate n {X — 0.8 in Fig. [5a|) . 
the population density of the intermediate species C in the 
(quasi-)steady state is reduced so much (cg — j ^ 0.01 in 
the figure) that the strong stochastic population oscilla- 
tions cause the evolution trajectories to touch the absorb- 
ing phase boundary; thus, the simulation eventually ter- 
minates at one of the absorbing states on the lines (a, 0, 6) 
which are not apparent in the mean-field analysis. 



Furthermore, as shown in Fig. I3b[ the quantitive re- 
lationship between the population density Cs for the in- 
termediate species C in the active coexistence state and 
the predation rate A numerically agrees remarkably well 
with the mean-field prediction Cs = fJ-fX. In Fig. [3cl we 
plot the temporal evolution of the three species' popu- 
lation densities for runs with predation rate A — 0.1; 
the population densities in the (quasi-)steady state (av- 
eraged over 50 runs at t = 3000 mcs) are found to be 
(0.4 ± 0.01, 0.1 ± 0.003, 0.4 ± 0.01), again fully consistent 
with the predictions of the mean-field approximation. We 
also numerically investigated the infiuence of quenched 
spatial disorder on the dynamical evolution of the sys- 
tem. However, in stark contrast with the robust fitness 
enhancement in the two-species LV model, in the hierar- 
chical three-species model spatial variability in the pre- 
dation rate A does not generate any noticeable effect on 
the dynamical evolution of the system. These observations 
are further elucidated by the spatial distribution snap- 
shot on a two-dimensional lattice depicted in Fig. I3dl All 
particles appear randomly scattered on the lattice in the 
(quasi-)steady state, and no correlated spatio-temporal 
structures such as either clusters or spirals become mani- 
fest. Consequently quenched spatial disorder cannot ben- 
eficially affect any correlations as in the LV model. 



Our Monte Carlo simulations demonstrate that the 
spatial three-species hierarchical food chain system, with 
both nearest-neighbor particle pair exchange and hopping 
processes present, is quantitatively well-described by the 
mean-field approximation. In the following subsection, we 
thus explore the respective role of particle pair exchange 
vs. hopping processes. 
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Fig. 3: ( Color online.) Observables for a stochastic three-species hierarchical food chain system in two dimensions 
with N — 256 x 256 sites, reaction rates fi — 0.01, cr — 0.1, and A as indicated; the simulations started with equal 
initial densities a(0) = b{0) — c(0) = 1/3 in (a) and (b), and a(0) — 6(0) = c(0) = 0.25 in (c) and (d); the data were 
averaged over 50 simulation runs; color coding in (d): blue/dark gray: A, red/gray: B, black: C, white: empty. 



4.2 Particle pair exchange processes and pure hopping 
processes 

In the above subsection, we have demonstrated that pop- 
ulation spreading through both nearest-neighbor particle 
pair exchange and hopping processes can well mix the 
two-dimensional spatially extended three-species hierar- 
chical food chain system. In order to further explore the 
origin for the appearance of mean-field behavior, we ran 
simulations wherein either the particle hopping or pair 
exchange processes were removed; Figs. Ha] and |4b] show 
snapshots of the resulting (quasi-)stationary particle dis- 
tributions. In Fig. l4al where only nearest-neighbor particle 
pair exchange was allowed, it is seen that the individuals 
for each species are uniformly distributed on the lattice, 
and no particular spatial patterns emerge in the system. 
We depict the associated temporal evolution (up to 6000 
mcs) of the three population densities for a smaller sys- 
tem {N = 100 X 100) in Fig. gel The population densi- 
ties (0.39 ± 0.01, 0.10 ± 0.01, 0.40 ± 0.02) in the (quasi-) 



steady state (averaged over 50 runs at t = 6000 mcs) are 
again quantitatively consistent with the mean-field pre- 
dictions. Also, as one would expect, we have not observed 
any remarkable influence of quenched spatial disorder on 
the evolutionary dynamics in these systems. 

Next, we investigate systems wherein particle spread- 
ing occurs solely through nearest-neighbor hopping, see 
Figs. |4b] and I4dl In Fig. I4bl one notices distinct species 
clusters, as also observed in two-species LV systems; how- 
ever, it is worthwhile noting that interactions between 
predators A and prey B do not directly take place at 
their respective cluster boundaries, but can only be indi- 
rectly realized via interspersed clusters of the intermediate 
species C (visible as black patches in Fig. I4bp . The re- 
quired presence of the intermediate species clusters signif- 
icantly coarse-grains the indirect reactions between preda- 
tors A and prey B, and effectively averages over spatially 
correlated structures. Again, we find that quenched spa- 
tial disorder does not noticeably affect the evolutionary 
dynamics of the system. However, as shown in Fig. l4d[ the 
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(a) Particle pair exchange processes 



(b) Nearest-neighbor hopping 
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Fig. 4: (Color online.) (a) and (b): Snapshots of the spatial particle distribution a.t t — 6000 mcs for a stochastic 
three-species hierarchical food chain system with N = 256 x 256 sites, reaction rates A = 0.1, fj, = 0.01, and a = 0.1, 
starting with equal initial densities a(0) = b{0) = c(0) = 0.25. Blue/dark gray: A, red/gray: B, black: C, white: empty, 
(c) and (d): Temporal evolution of population densities for simulation runs with (a) only nearest-neighbor particle 
exchange, and (b) only nearest-neighbor hopping, with N = 100 x 100 sites, averaged over 50 simulation runs. 



population densities (0.29 ± 0.01, 0.14 ± 0.01, 0.47 ± 0.02) 
in the (quasi-) steady state (averaged over 50 runs at 
t — 6000 mcs) do not coincide with the predictions from 
the mean- field rate equations, i.e., nearest- neighbor hop- 
ping processes alone are not strong enough to effectively 
mix the two-dimensional three-species hierarchical food 
chain system. We have furthermore checked that these 
numerical observations, here obtained for A = 0.1, about 
the effect of either particle pair exchange processes or pure 
hopping processes also apply to systems with much smaller 
predation rate A — 0.02. 

In summary, we have demonstrated that the hierarchi- 
cal three-species food chain model with site restrictions 
and with particle spreading through pair exchange pro- 
cesses effectively washes out the formation of any spatial 
structures, hence becomes well mixed, and behaves as in 
the mean-field predictions. As is also observed in the two- 
species LV system with site restrictions, the hierarchical 
food chain model also possesses a critical threshold for the 
predation rate, where a phase transition from an inactive 



absorbing state to an active coexistence state occurs. How- 
ever, if particle spreading happens solely through nearest- 
neighbor hopping, distinct species clusters emerge. Yet 
since intermediate species C clusters are necessary to indi- 
rectly promote predation interactions between the A and 
B species, the effective reaction boundaries between the 
predators and prey become blurred, such that any spatial 
variability still cannot remarkably infiuence the dynami- 
cal evolution of these food chain systems. Consequently, 
at least for the parameter values and system sizes stud- 
ied here, the ensuing hierarchical three-species food chain 
model does not closely resemble the features of the corre- 
sponding two-species spatial stochastic LV system that is 
governed by strong fluctuations and marked correlations 
in the active coexistence state. 
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5 Conclusion 

In this paper, we studied the connection between two 
three-species predator-prey systems and the well-understood 
two-species Lotka-Volterra model. First, we explored the 
evolutionary behavior of the two minority species in the 
"corner" three-species cyclic rock-paper-scissors (RPS) 
model with strongly asymmetric rates. We analytically 
demonstrate that in the mean-field limit the evolution- 
ary dynamics of both minority species in the "corner" 
RPS model can be well approximated by the two-species 
LV model. Employing Monte Carlo simulation for a two- 
dimensional spatially extended "corner" RPS model, we 
found the population densities for those two minority species 
in the (quasi-)steady state to coincide with both the an- 
alytic predictions from mean-field approximation in the 
resulting two-species LV model and its associated numeri- 
cal simulation results. Furthermore, introducing quenched 
spatial disorder into the predation rate A in the "corner" 
RPS model, we observe that the influence of quenched spa- 
tial disorder on the evolution of system cannot be ignored, 
in contrast to RPS systems far away from the "corner" 
of configuration space, wherein quenched spatial disorder 
has only minor effects on the evolution of the system [52., : 
The minority species' fitness in the "corner" RPS system 
is markedly enhanced due to the spatial variability in the 
predation rates, as is indicated by larger population den- 
sities for both minority species in the (quasi-)steady state, 
and accompanied by shorter relaxation times and stronger 
localization of the clusters for those two minority species. 
That is, the two minority species in the "corner" RPS 
model behave just like the predator and prey in the spa- 
tially extended two-species LV system. 

We also investigated hierarchical three-species "food 
chain" systems, in which an additional intermediate speci- 
es is inserted between the predators and prey in the classic 
two-species LV model. As observed in the spatial version 
of the two-species LV model with site occupancy restric- 
tions, in the similarly restricted three-species "food chain" 
model with finite total carrying capacity, there exists a 
critical predation rate threshold, which in the thermody- 
namic limit represents a non-equilibrium phase transition 
from an inactive absorbing to an active coexistence state. 
Our simulation results show that nearest-neighbor parti- 
cle pair exchange processes are sufficient to wash out the 
formation of species clusters and consequently generate a 
well-mixed spatial system, whose features can be quantita- 
tively well approximated by the mean-field rate equations, 
and wherein the influence of quenched spatial disorder 
may safely be ignored. Moreover, if we allow only nearest- 
neighbor hopping processes, and eliminate any pair ex- 
changes, distinct species clusters appear in the coexistence 
state. However, due to the interruption by intermediate 
species clusters, the predation processes between preda- 
tor and prey are effectively coarse-grained, and fluctuation 
and correlation effects on the resulting two-species dynam- 
ics suppressed. Consequently, quenched spatial disorder 
still cannot enhance the fitness of either species. That is, 
as opposed to the "corner" rock-paper-scissors model with 
strongly asymmetric rates, the spatially extended hierar- 



chical three-species "food chain" system actually does not 
behave like the two-species LV model. 
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